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1. INTRODUCTION 

In solving the transonic fhll potential equations by using finite 
difference approximations on a transformed coordinate system, the solution 
is affected by the following four kinds of numerical errors: 

(1) numerical error in evaluating the metrics of the coordinates, 

(2) error due to discretization of the full potential equation, 

(3) error due to artificial viscosity for the supersonic domain, and 

(4) error due to approximating the shock behavior on the grids 
with finite spacing. 

When a numerical grid generation method is used, the matrics of the 
transformation are calculated by finite difference approximations using 
the coordinates of the grids, so they are subject to errors, which in turn 
affect the accuracy of the finite difference approximation of the flow 
equation on the transformed coordinates. The error due to discretization 
is attributed to the truncation error of the finite difference 
approximation of the flow equation and becomes larger as the partial 
derivatives of the flow potential distribution become larger. The 
artificial viscosity is used to stabilize the iterative convergence of the 
transonic full potential equations, but introduces errors in the solution. 
Although the error in the numerical solution and the truncation error are 
related to each other in a very complicated manner, the error of the 
solution decreases when the truncation errors of the difference 
approximation in both the metric calculation and the fl"w solution 
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decrease. Naturally it is desirable to increase the number of grids in 
the areas of the flow domain where the truncation error is likely to be 
high. 

Recent studies [1-31 on the use of adaptive grids for transonic full 
potential flow computations indicate that computing efficiency in terms of 
accuracy and cost is substantially improved by using adaptive grids. In 
those studies, the computations proceed in the following four steps: (1) 
generation of initial grids, (2) solution of the full potential equation 
on the initial grids (initial flow solution), (3) generation of adaptive 
grids using the information from the initial flow solution, and (4) flow 
solution on the adpative grids (adaptive grid flow solution). In 
generating adaptive grids, the locations of the supersonic bubble, the 
location of shock, and the area of high flow acceralation are identified 
from the initial flow solution. 

In order to further develop the adaptive grid procedure, however, two 
major questions should be answered. First, what information from the 
solution on the initial grids should be used in generating the adaptive 
grids, and second, how to allocate the total number oi grids available in 
the entire flow domain in such a way that the overall accuracy of the 
solution becomes optimized. Although these problems were investigated in 
recent studies [4,5], no clear answers have been obtained in conjunction 
with the full potential flow calculations. 

The objectives of the present study are (1) to investigate the 
feasibility of using an estimate of the truncation error of the difference 
equation as a measure for the quality of the grids, (2) to investigate 
what information from the initial solution should be used in generating 
adaptive grids, and (3) the effect of irregularity of grid spacings. 

For the s^ke of simplicity, our effort in this study is 
focused on the grid distribution along the airfoil surface. Furthermore, 
the numerical investigation is performed by solving the full potential 
equation on the Cartesian coordinates for a symmetric parabolic arc 
airfoil. Use of coordinate transformation is avoided so that the study 
can be done without being affected by the numerical error of the metrics 
of the transformation. Therefore, the airfoil boundary condition in the 
numerical study is represented along a straight line as a thin airfoil 
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perturbation. 

2. TRUNCATION ERROR ALONG THE AIRFOIL SURFACE 

In this section, the effect of grid spacing on the truncation error 
of a finite difference approximation for the two-dimensional full 
potential equation on the Cartesian coordinates is analytically studied. 
The full potential equation is written as 

(R(F x»x + <« p y »y s 0 (D 

where F is the potential fuction and R is the fluid density. We proceed 
the derivation of the truncation error expressions for the following 
cases: (A) R=1 (flow is assumed to be incompressible) and constant grid 
spacing, (B) R=1 (incompressible flow) and variable grid spacing, (C) 
compressible subsonic flow with variable grid spacing (R is variable) , 
and (D) supersonic flow with variable grid spacing. In the remainder of 
this report, only the truncation error of the difference approximation for 
the first term of Eq.(1) along the airfoil surface will be considered. 

(A) Truncation Error for Incompressible Flow on Equally 
Spaced Grids 

With R=1 and equally spaced grids on the x-coordinate, the difference 
approximation for the first term of Eq.(1) is written as 

«W-. (!) 

By applying the Taylor expansion to F^ j and j, Eq.(2) 
becomes 

fa/6x P - s F xx + (1/12)h 2 F iv ♦ (1/360)hV l ♦ ... (3) 

where h is the grid spacing and F n is the n-th partial derivative of F 
with respect to x. Equation (3) shows that the truncation error is 
represented by 

TE(x) = d?F/6x 2 - F XJ( * ( 1 / 12 )h 2 F iv ' (M) 
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which is second-order in h. 

(B) Truncation Error for Incompressible Flow on Nonuniform Grids 

Using the three grids at x-g, x, and x+h, respectively, the 
three-point second order difference approximation for the first term of 
Eq.O) is written as 


, , [F(x+h)-F(x)Vh - fF(x)-F(x-g)1/g 

fi“F/Sx z = 

(g+h)/2 


(5) 


Using the Taylor series expansion, the truncation error is written as 
TE(x) * 5 2 F/55? - F xx 

<(g/h) 3 ♦ 1 
g/h ♦ 1 


* 0/3)(h-g)F «• ♦ <h 2 /l2) — F iv ♦ ... (6) 


The above equation shows that, if gjih, the truncation error is 
proportional to (h-g), indicating that a rapid change of grid spacing 
causes a significant amount of truncation error. 

(C) Truncation Error for Subronic Compressible FLow on 
Nonuniform Grids 


The difference approximation for the first term of Eq.O) on 
nonuniforaly spaced grids may be written as 


<g-h)/2 r Ri+1/2 



R i-1/2 



(?) 


where g and h are the grid spacings defined in (B). The truncation error 
of Eq.(7) may be found by introducing Taylor series expansion of R^ ^ 
and ^ about x^ into Eq.(7) and subtracting the first term of 
Eq.O). After a manipulation of the terms, the lowest terms in the 
truncation error are found as 
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TE(x) = 


♦ (h-g) [(1/l65R i »’F i ’ ♦ (1/iOR i 'F i t * ♦ (1/3)R 1 F i ”*3 


g^+h^ 


g 2 +h 2 


gW 


♦ [ R" *F* ♦ R' ‘F’ • R'F" ’ 

9S(g+h) 8 3(g+h) 

g^+h^ 

+ R F" "] 

12(g+h) 


( 8 ) 


where R', R*', and R ,M are partial derivatives of R with respect to x. 
The first term of Eq.(8) is proprtional to (h-g) and vanishes vften g=h. 

(D) Truncation Error for Supersonic Flow on Nonuniform Grids 


The full potential equation for the supersonic domain is written as 

(K< F X ) X ♦ (R(F) > s 0 (q) 

where R is the upwind approximation for R in the supersonic domain. The 
upwind shifting of density is adopted to stabilize the iterative scheme 
for Eq.(9). Upwind density is given as 

*i-1/2 = wR i-3/2 * (1 ’ w)R i-l/2 (1C> 

where w is a upwind-shiftirg parameter. 

By writing 


*i-3/2 s R i— 1/2 “ AxR ’i-1/2 

where 

AX = x i-1/2 - x i _ 3/2 

Eq.(10) becomes 

^i-1/2 = R i-1/2 ” w<lxR 'i-1/2 
Then, the first term of Eq.(9) becomes 


( 11 ) 


ff (F x))„ * (R(F X » X - wax(R'(F x )) x 
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Since the first term of Eq.(12) is identical to that analyzed in (C) , the 
truncation error of the first term is given by Eq.(B). The second term on 
the right side of Eq.(12) is the truncation error caused by upwind 
shifting of fluid density. The total truncation error of the supersonic 
domain is therefore written as 

TE(x> * TE q (x) ♦ TE UVJ (x> (13) 

where TE Q (x) is the truncation error that occurs vhen no upwind shifting 
of density is used and given by Eq.(B), and TE uy (x) is equal to the 
second term of Eq.(12): 

^u u (x) s -WAX(R'(F X )) X (1U) 

3 . NUMERICAL EVALUATION OF TRUNCATION ERROR AND OTHER 
QUANTITIES OF THE aCM SOLUTION ON THE INITIAL GRIDS 

Numerical evaluation of the truncation error formulas derived above 
is of interest because they may be used as means to evaluate the quality 
of the grids used and as grid spacing control (GX) function. It should 
be noticed, however, that the formulas derived above are expressed in 
terms of the exact flow solution. In the practice of adaptive grid 
generation and the flow solution on them, the only way of the numerical 
evaluation is to use, in the truncation error equations, an numerical (or 
approximate ) solution on the initial grids in place of the exact 
solution. Therefore, the trncation error estimated by using the finite 
difference solution is subject to the error of the solution which are 
caused by the truncation errors of on the initial grids. Furthermore, the 
estimated truncation error is influenced by the nonlinear nature of the 
full potential equation where R is a function of the solution. (It should 
be noticed that, in the previous section, R is assumed to be known 
exactly.) 

In spite of this limitation we attempt to estimate the truncation 
error of the ntmerical solution of the flow solution. Instead of using 
the truncation equations derived above, the truncation is evaluated 
directly by 
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TE(x) s 6RSF/6X 2 - 3R3F/3x 2 (15) 

where the second term is evaluated by the 5-point difference approximation 
which is fourth-order-accurate. 

The full potential equation was first solved on the Initial grids, 
plotted in Fig. 1, in which 41 grids are uniformly spaced along the wing. 
The airfoil was assumed to be a parabolic arc airfoil with 10J thickness 
ratio. In order to avoid the use of a coordinate transformation, the 
airfoil was represented by the boundary condition on the bottom straight 
line written as 

3F/3y s f(x)M (16) 

where f(x) is the surface profile of the airfoil, and M is the Mach 
number. Equation (15) is a crude approximation of the boundary condition, 
but sufficient as a computational model fpr the present study since all 
the calculations hereafter are performed by using the same model. The 
difference equations for the full potential equations were iteratively 
solved by using the AF1 scheme [6]. 

The following quantities of the initial flow solution, including the 
truncation error, all along the wing surface were calculated and plotted 
in Fig. 2. 

a) CP (pressure coefficient) distribution 

b) An estimate for the truncation error of the difference 
approximation 

c) Second derivative of the flow potential distribution 

d) Truncation error due to upwind shifting of fluid density, 

!TE uw (x)I 

e) The space derivative of the fluid density along the airfoil 
surface 

Figure 2-A is the CP distribution. Figure 2-3 is the truncation error 
JTE(x)i along thj airfoil surface. ',TE(x)5 jumps up at the grid next to 
the leading edge, but decreases quickly as x i z increased until xr0.25 is 
reached, and then increases because of the effect of upwind shifting of 
the flow density. The peaks around xs0.7 are due to the shock. After 
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reaching a bottom at around xs0.8, ',TE' increases toward the trailing edge 
becuase of the increasing negative gradient of the airfoil. Figure 2-C 
shows ',F W „! along the airfoil surface. jF ! has a great similarity 
to JTE! except in approximately 0.25<x<0.? where the supersonic bubble and 
shock are located. Figure 2-D shows the truncation error due to the 
upwind snifting of the density. Figure 2-E is 13R/3x! along the airfoil, 
which shows its high sensitivity to the location of shock. 

Very small values of |TE(x)| of Fig. 2-9 in the range 0.1<x<0.25 and 
around x=0.8 seem to indicate tf at in these areas much greater grid 
spacing than in the initial grids may be used. The two peaks of ITE(x)', 
adjacent to the leading and trailing edge, which are due to the highest 
curvature of the airfoil at the edges, are indicative of the need of finer 
grids toward the edges. The two highest peaks at about xs9.7 due to the 
shock can not be interpreted simply as the truncation error because of the 
strong nonliiear nature of the full potential equation across the shock, 
but should be rather interpreted as an indicator of the shock location 
which demands finer grids. 

One might ask if the adaptive grids that will yield a levelized 
JTE(x) ' can be considered as an optimized adaptive grids (or at least as 
better grids than the initial grids). A part of the answer will be sought 
in the next section. 


4. EFFECT OF DIFFERENT GRID SPACING CONTROL FUNCTIONS ON 
ADAPTIVE GRIDS AND FLOW SOLUTIONS 

Effects of using different quantities obtained from the initial flow 
solution as the GDC function, with or without modifications, on the 
adaptive grids generated and the flow solution are studied in thi3 
section. The procedure of generating adaptive grids along the airfoil 
surface is first described. 

Assuming that only the grids along the airfoil surface are 
regenerated by the adaptive grid generation procedure, two functions are 
defined as 

S(x): grid density control (GDC> function 



<*#»=. «*»$>* , , 
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iKx): grid density function 

where x is the coordinate along the airfoil surface. The GDC function is 
a function that carries such information of the initial flow solution as 
the locations of the supersonic bubble, the shock, and areas of high 
acceleration and deceleration. The grid density function is a function 
that specifies the density of grids as a function of x, and is inversely 
proportional to the grid spacing. 

The two functions defined in the above paragraph are related by 
d 2 

S”’y ( X ) ♦ 0 !j>(x) = S(x) ( 17 ) 

dx 2 

i|>’(0) = 0'(1) = 0 

where 0 is a control parameter, x=0 and xsl are the leading edge and 
trailing edge of the airfoil, respectively. The above equation is solved 
by finite difference approximations using the tridiagonal solution scheme 
once S(x) is specified. Then, the coordinates of the grids along the 
airfoil surface are obtained by 

x k s T^Ck/K), ksO,1,2,...K 08) 

where K is the total number of grid intervals on the airfoil surface, and 
k is a grid number counted from the leading edge. In the above equation, 
T(x) is given by 


x 

T(x) = c f «(s)ds 09) 

J 0 

where e is a normalizing constant to let T(1)=1. 

Equation 06) may be considered as a tranformation of S(x) into <l<x^. 
The purpose of this transformation is to smooth out any abrupt change in 
S(x) in order to avoid a sudden change of spacing of adaptive grids. The 
control parameter a determines the intensity of the smoothening effect. 

As o is increased, * becomes more and more similar to S(x>, thus the 
3moothening effect disappears. On the other hand, a very small value will 
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information carried by S(x). 

In the remainder of this section, the results of five case studies 
using different GDC functions are discussed. 

Case 1 

The adaptive grids were generated by using STEu)! as the GDC (grid 
density control) function shown in Fig. 2-9, on which the flow equation 
was solved. In Figure 3, the CP distribution on the adaptive grids are 
compared with that on the initial grids. The adaptive grids generated 
along the wing are shown in Fig. 3 by tick marks along the horizontal 
coordinate line at the bottom of the figure. 

In the adaptive grids thus generated, grids are strongly clustered 
toward the leading and trailing edges as well as around the shock. On the 
other hand, the grids (in 0.1<x<0.25) that are between the leading edge 
and the supersonic region and the grids (in 0.B<x<0.9) between the shock 
and the trailing edge became coarse. The strong clustering in the three 
areas are resulted by the peaks of !TE(x)! of the flow solution on the 
initial grids near the leading and trailing edges as well as at the shock. 
At the same time, the coarse grids in the other two areas are due to very 
small values of |TE(x)l as seen from Fig.2-B (see the line witn lozenge 
marks) . 

A comparison of the two CP distributions in Fig. 3 suggests several 
effects of the grid spacings • n the flow solution. Coarser grids in the 
subsonic region after the leading edge causes more expansion of the fluid 
than on finer grids (CP distribution on the adaptive grids in 9.1<x<0.25 
are higher than that on the initial grids, where the adaptive grids are 
coarser than the initial grids). The CP distribution at the shock on the 
adaptive grids are much sharper than that on ‘■he initial grids. The 
expansion of the fluid before the shock in the range of 0.5<x<0.6 on the 
adaptive grids are suppressed compared with that on the initial grids. 

This is attributed to the coarser grids in this a**ea of the adaptive grids 
than those on the initial grids. The coarser grids in the supersonic 
domain results in using a higher fluid density in the difference equation 
because of the upwing shifting of the density in the supersonic region. 
However, fine grids just before the shock caused an opposite effect of the 
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upwind shifting that results in a faster expansion of the fluid. 

Case 2 

}TE(x) | is used as the GX function just as in Case 1 except it is 
modified so as to decrease the strength of clustering and to increase the 
grid density in the coarse grid area. The GDC function is defined as 
follows: 

S(x) a !TE(x> * if 0.5<!TE(x)!<0.0; 

S 0.5 if 'TE(x) ! >0.5 

a 0.01 if JTE(x) { <0.01 (20) 

The CP distribution on the adaptive grids is compared with that of the 
initial solution in Fig. H-A. Although expansion and contraction of Rrids 
became much milder than in the pr : us case, it is still seen that (1) a 
rapid expansion of grid spacing around 0.1 to 0.25 causes an increase in 
CP on the graph, and (2) a rapid contraction of grid spacing toward the 
shock causes a lower CP on the graph. Figure 4-B shows a comparison of 
STE(x)! between the initial solution and the adaptive grid solution. It 
is noticed from this figure that JTE(x)! in the range of 0.1<x<0.25 for 
the adaptibe grid flow solution is still very small despite the 
discrepancy in CP. This indicates that the |TE(x)| is insensitive to the 
change of physically important quantities of the flow solution. 

The results of both Case 1 and the present use 2 show that estimated 
truncation error distribution is a poor indicator of the error of the 
solution, particularly for the error of the CP distribution. 

Case 3 

The adaptive grids are generated using !F xx (x'i shown in Fig. 2-C 
as the GX Junction. Since !F XX I has sharp peaks near the leading and 
trailing edges, the grids shown in Fig. 5 are strongly clustered toward 
the leading and trailing edges. The IF xx (x)| has a pe3k at the shock, 
which causes a clustering at the shock. The grid density at the shock 
becomes even coarser than that on the initial grids ct the shock because 
of much stronger clustrring effects at the leading and trailing edges. 

The CP distribution on the adaptive grids shown in Fig. 5 loorcs poor. 
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The same effect of coarseness on the CP distributions as seen before are 
obverved. A rapid expansion of grid spacing in the subsonic domain causes 
an upward shifting the the CP distribution, while coarse grid in the 
supersonic domain causes a downward shifting of CP and smearing of the 
shock. 

The use of |F xx (x>{ results in a poor accuracy as a whole because 
the strong clustering toward the leading and trailing edges causes a lower 
grid density in all other areas. 

Case 4 

i3R/3x| is used as the GDC function. The CP distribution on this 
choice of the GDC function and the adaptive grids generated are shown in 
Fig. 5. With this choice, the adaptive grids are moderately clustered 
toward the leading and trailing edges, while the grids around the shock 
are strongly clustered. No particular clustering is seen in the 
supersonic bubble. Ccaser grids In both subsonic and supersonic 
expansion region causes the same effects as seen in the previous cases. 

Case 5 

The GDC function is set as 

$(x) s M ? (x) (211 

where M is the Mach number, except modified at the leading and trailing 
edges, and the shock as 

S(x) = 2M ? (x) (22) 

The purpose of the modification by Eq.(22) is to enhance the clustering 
effect at the shock. Figures 7-A, 8-A, and 9-A are the CP distributions 
on the adaptive grids using 41, 31 and 21 grids along the airfoil surface, 
respectively. They are compared to the CP distribution on the initial 
grids. Each of Fig. 7-B, 8-B and 9-B shows the ITE(x)! distributions 
corresponding to the two CP distributions in Figs. 7-A, R-A and 9-A, 
respectively. 
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The grids generated with the present choice of the GDC function 
yields smooth changes of grid spacing. Each of Figs. 7-A, 8-A, and 9-A 
shows that the shock becomes better represented on the adaptive grids 
than on the equally spaced initial grids, while coarser grids upwind of 
the shock little affect the CP distribution. The disagreement of the CP 
distribution downwind of the shock between the initial flow solution and 
the adaptive grid flow solution is attributed to the change of the CP 
distribution at the shock. 

A comparison of the CP distributions of Fig. 7-A to those in Fig. 

9-A is instructive. Even though the grids in Fig. 9-A are coarser than 
those on Fig. 7-A, the CP distributions upwind the shock agree well. 
Although the coarser grids in Fig. 9-A causes more expansion upwind of 
the shock than in Fig. 7-A, the difference is smaller than the 
discrepancies caused by a rapid change of grid spacing as seen in the 
previous cases. 

5. EFFECT OF GRID SPACING IRREGULARITY ON THE FLOW SOLUTION 

During the case studies in the previous section, it was identified 
that a sudden change of grid spacing may cause a kink in the CP 
distribution. This is considered as the effect of the first order 
truncation error in Eq.(S) due to nonuniforrity of grid spacing. 

In order to study the effect of irregularity of the grids 
quantitatively, the flow solution was repeated with equally spaced grids 
except th3t only one of the grids along the airfoil surface is moved 
backward or forward at a time. In all calculations for this study, only 
21 grids are used along the airfoil. 

Figure 10 shows two C? distributions, where the third grid counted 
from the leading edge is moved half way forward or backward to the next 
grid as marked by lozenge and circular marks at the bottom of Fig. 10. 
This figure shows th3t moving a grid backward (see lozenge marks) causes 
an upward deviation of the CP curve on the moved grid while the CP on the 
adjacent grid upwind is deviated downward. 

Figure 11 shows similar results of moving the 10th grid, which is 
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located in the supersonic flow domain. The deviations caused by the move 
show almost the same trend as seen in Fig. 10, namely the forward shift of 
a grid causes a a upward shift of CP on the graph at that grid, while the 
backward shift causes a downward shift of CP on the graph. 

Figure 12 shows the results for the 18th grid moved. The deviation 
of the CP distribution caused by the move is opposite to those seen in 
Figs. 10 and 11, suggesting that the direction of the CP deviation depends 
on Aether the grid is in. a expanding flow or in a contracting flow. 


6. CONCLUSIONS 

This paper described the results of the study on (i) the analytics? 
expressions for the truncation error of the transonic full potential 
equation, (ii) nir.-arical study of the estimated truncation error, (iii) 
adaptive grid generation using various different quantities obtained from 
the initial flow solution, and (iv) effect of irregularity of grid spacing 
on the flow solution. 

The estimated truncation error does not seem to represent the error 
of the flow solution adequately, thus is not suitable as a grid density 
control (GDC) function for adaptive grid generation. Among several 

P 

different quantities taken from the initial flow solution, the M 
distribution is found to be most useful with a modification to enhance the 
clustering at the shock. 

One important fact learned quantitatively through the numerical study 
is that a rapid change of grid spacing is damaging to the accuracy of the 
flow solution. This means that, while the objective of adaptive grids is 
to dispatch a larger number of grids to where the error is high, the grid 
density should be only gradually changed from the low error area to the 
high error area. The need of gradual change of grid density justifies the 
use of the diffusion process for the GX function proposed by the author. 

The practice of r mathematical optimization of grid distribution for 
the full potential euq< »ion does not seem to feasible, because of the 
complex nature of the relation among an error estimate, grid spacing and 
the real error of the flow solution. For further research, it is rather 
recommended to develop empirical relations between the CP distribution and 


lU 



ORIGINAL PAGE 13 
OF POOR QUALITY 


the required grid spacing distribution as well as the acceptable rates of 
grid spacing change. 
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Figure 2-E Distribution of|^P. / e>x| 
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Figure 4-A Pressure Coefficient., Case 2 
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Figure 5 Pressure Coefficient., Case 3 
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